{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kim2014 identification statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "\n",
    "src_dir = os.path.abspath('../src')\n",
    "if src_dir not in sys.path:\n",
    "    sys.path.append(src_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import functools\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as mticker\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import tqdm\n",
    "\n",
    "from ann_solo import reader, spectrum, utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot styling\n",
    "plt.style.use(['seaborn-white', 'seaborn-paper'])\n",
    "plt.rc('font', family='serif')\n",
    "sns.set_palette('Set1')\n",
    "sns.set_context('paper', font_scale=1.)    # single-column figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "tqdm.tqdm = tqdm.tqdm_notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_time_from_log(filename):\n",
    "    with open(filename, 'r') as f_in:\n",
    "        for line in f_in:\n",
    "            if line.startswith('real'):\n",
    "                # Wall clock time.\n",
    "                realtime = line.split()[1]\n",
    "                minutes = int(realtime[:realtime.find('m')])\n",
    "                seconds = float(realtime[realtime.find('m') + 1:\n",
    "                                         realtime.rfind('s')])\n",
    "                realtime = minutes * 60 + seconds\n",
    "                \n",
    "                return realtime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_fdr = 0.01\n",
    "tol_mass = 0.1\n",
    "tol_mode = 'Da'\n",
    "min_group_size = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "id_dir = '../../data/processed/kim2014/gpu'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "edd2317207c14fb6b28a0950df419121",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='Files processed', max=4207, style=ProgressStyle(description_w…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "ssms = []\n",
    "num_ssms, num_peptides, runtimes = [], [], []\n",
    "for filename in tqdm.tqdm(os.listdir(id_dir), desc='Files processed',\n",
    "                          unit='files'):\n",
    "    filename_full = os.path.join(id_dir, filename)\n",
    "    filename_base, ext = os.path.splitext(filename)\n",
    "    if ext == '.log':\n",
    "        runtimes.append((os.path.splitext(filename)[0],\n",
    "                         extract_time_from_log(filename_full)))\n",
    "    elif ext == '.mztab':\n",
    "        file_ssms = reader.read_mztab_ssms(filename_full)\n",
    "        ssms.append(file_ssms)\n",
    "        file_peptides = (file_ssms['sequence'].str.replace(r'n?\\[\\d+\\]', '')\n",
    "                         .unique())\n",
    "        num_ssms.append((filename_base, len(file_ssms)))\n",
    "        num_peptides.append((filename_base, len(file_peptides)))\n",
    "        \n",
    "ssms = pd.concat(ssms, ignore_index=True)\n",
    "\n",
    "num_ssms_df = pd.DataFrame.from_records(num_ssms,\n",
    "                                        columns=['filename', 'ssms'])\n",
    "num_peptides_df = pd.DataFrame.from_records(num_peptides,\n",
    "                                            columns=['filename', 'peptides'])\n",
    "time_df = pd.DataFrame.from_records(runtimes, columns=['filename', 'time'])\n",
    "\n",
    "stats = functools.reduce(\n",
    "    lambda left, right: pd.merge(left, right, on='filename'),\n",
    "    [num_ssms_df, num_peptides_df, time_df])\n",
    "stats = stats.sort_values('filename').reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of identified spectra: 14,032,494\n",
      "Total search time: 281 h (11.7 days)\n",
      "Search time per file: 8.0 ± 3.1 min\n"
     ]
    }
   ],
   "source": [
    "print(f'Number of identified spectra: {stats[\"ssms\"].sum():,.0f}\\n'\n",
    "      f'Total search time: {stats[\"time\"].sum() / 3600:.0f} h '\n",
    "      f'({stats[\"time\"].sum() / 86400:.1f} days)\\n'\n",
    "      f'Search time per file: {stats[\"time\"].mean() / 60:,.1f} '\n",
    "      f'± {stats[\"time\"].std() / 60:,.1f} min')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_ssm_unmod = num_ssm_total = 0\n",
    "for filename in os.listdir(id_dir):\n",
    "    filename_full = os.path.join(id_dir, filename)\n",
    "    filename_base, ext = os.path.splitext(filename)\n",
    "    if ext == '.log':\n",
    "        with open(filename_full, 'r') as f_in:\n",
    "            for line in f_in:\n",
    "                if line.endswith('spectra identified after the standard search\\n'):\n",
    "                    start = line.find(' : ') + 3\n",
    "                    end = line.find(' ', start)\n",
    "                    num_ssm_unmod += int(line[start:end])\n",
    "                if line.endswith('spectra identified after the open search\\n'):\n",
    "                    start = line.find(' : ') + 3\n",
    "                    end = line.find(' ', start)\n",
    "                    num_ssm_total += int(line[start:end])\n",
    "num_ssm_mod = num_ssm_total - num_ssm_unmod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of identified spectra: 14,032,494\n",
      "Number of unmodified identified spectra: 9,760,497\n",
      "Number of modified identified spectra: 4,271,997\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(f'Total number of identified spectra: {num_ssm_total:,.0f}\\n'\n",
    "      f'Number of unmodified identified spectra: {num_ssm_unmod:,.0f}\\n'\n",
    "      f'Number of modified identified spectra: {num_ssm_mod:,.0f}\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_mass_groups(ssms, tol_mass, tol_mode):\n",
    "    ssms_remaining = ssms.sort_values('search_engine_score[1]',\n",
    "                                      ascending=False)\n",
    "    ssms_remaining['mass_diff'] = ((ssms_remaining['exp_mass_to_charge'] -\n",
    "                                    ssms_remaining['calc_mass_to_charge']) *\n",
    "                                   ssms_remaining['charge'])\n",
    "\n",
    "    # Start with the highest ranked SSM.\n",
    "    mass_groups = []\n",
    "    while ssms_remaining.size > 0:\n",
    "        # Find all remaining PSMs within the mass difference window.\n",
    "        mass_diff = ssms_remaining['mass_diff'].iloc[0]\n",
    "        if (tol_mass is None or tol_mode not in ('Da', 'ppm') or\n",
    "                min_group_size is None):\n",
    "            mask = np.full(len(ssms_remaining), True, dtype=bool)\n",
    "        elif tol_mode == 'Da':\n",
    "            mask = (np.fabs(ssms_remaining['mass_diff'] - mass_diff) <=\n",
    "                    tol_mass)\n",
    "        elif tol_mode == 'ppm':\n",
    "            mask = (np.fabs(ssms_remaining['mass_diff'] - mass_diffs) /\n",
    "                    ssms_remaining['exp_mass_to_charge'] * 10 ** 6\n",
    "                    <= tol_mass)\n",
    "        mass_groups.append(ssms_remaining[mask])\n",
    "        # Exclude the selected SSMs from further selections.\n",
    "        ssms_remaining = ssms_remaining[~mask]\n",
    "\n",
    "    mass_group_stats = []\n",
    "    for mass_group in mass_groups:\n",
    "        mass_group_stats.append((mass_group['mass_diff'].median(),\n",
    "                                 mass_group['mass_diff'].mean(),\n",
    "                                 len(mass_group)))\n",
    "    mass_group_stats = pd.DataFrame.from_records(\n",
    "        mass_group_stats, columns=['mass_diff_median', 'mass_diff_mean',\n",
    "                                   'num_ssms'])\n",
    "    return mass_group_stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "mass_groups = get_mass_groups(ssms, tol_mass, tol_mode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mass_diff_median</th>\n",
       "      <th>mass_diff_mean</th>\n",
       "      <th>num_ssms</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000646</td>\n",
       "      <td>0.001236</td>\n",
       "      <td>9882777</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>57.022315</td>\n",
       "      <td>57.022839</td>\n",
       "      <td>308387</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>27.995986</td>\n",
       "      <td>27.998766</td>\n",
       "      <td>246428</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>0.993853</td>\n",
       "      <td>0.987444</td>\n",
       "      <td>219006</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>15.995468</td>\n",
       "      <td>15.996442</td>\n",
       "      <td>211927</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>-0.986346</td>\n",
       "      <td>-0.999017</td>\n",
       "      <td>163269</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>14.015814</td>\n",
       "      <td>14.013037</td>\n",
       "      <td>133020</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>76</th>\n",
       "      <td>-17.025133</td>\n",
       "      <td>-17.019578</td>\n",
       "      <td>129687</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>54</th>\n",
       "      <td>-18.009568</td>\n",
       "      <td>-18.006636</td>\n",
       "      <td>111075</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.988133</td>\n",
       "      <td>1.989208</td>\n",
       "      <td>99286</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    mass_diff_median  mass_diff_mean  num_ssms\n",
       "0           0.000646        0.001236   9882777\n",
       "8          57.022315       57.022839    308387\n",
       "21         27.995986       27.998766    246428\n",
       "16          0.993853        0.987444    219006\n",
       "11         15.995468       15.996442    211927\n",
       "5          -0.986346       -0.999017    163269\n",
       "7          14.015814       14.013037    133020\n",
       "76        -17.025133      -17.019578    129687\n",
       "54        -18.009568      -18.006636    111075\n",
       "4           1.988133        1.989208     99286"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mass_groups.sort_values('num_ssms', ascending=False).head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqcAAADQCAYAAAA6XchUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X1cVHX+//8HoOAA0qBoKF6gpqiQmGuMbWymUu26u5X7sbZvkdW6rNbG5sUKa2t+unBTMK9RuthuWu56kSnaxa67a+a2GmBWJiQiLRQi4K4xg8iQg3B+f/hzPk6AoHExwPN+u3G7MWfe58zrnDPDPHmf8z7HwzAMAxERERERN+DZ1gWIiIiIiFykcCoiIiIibkPhVERERETchsKpiIiIiLgNhVMRERERcRsKpyIiIiLiNhRORURERMRtKJyKiIiIiNtQOBURERERt6FwKiIiIiJuQ+FURERERNyGwulVCgsLa+sSRERERDochVMRERERcRsKpyIiIiLiNhRORURERMRtKJyKiIiIiNtQOBURERERt6FwKiIiIiJuQ+FURESknfHw8GjrEkRaTJe2fPGSkhIWL17MqFGjOHLkCI899hiDBw8mOTmZoKAgrFYrCQkJeHl5kZOTw9atW+nWrRtjx44lJibGZVmGYbB8+XL8/f0pKSkhMTERk8lEcXExL774Ij179qRXr17cf//9AOzcuZOCggLKysqYPn06oaGhbbAFRERERORSbdpzWl1dzb333ssvf/lL7rzzTl5//XXS0tIIDQ1l5syZmEwm9u7dC8CiRYuYM2cOCQkJrFu3DofD4bKsAwcOcP78eWbMmEF4eDhbt24FYNmyZcTGxvLEE0+we/duTp8+TUVFBW+88QazZ88mLi6OpKSkVl93EREREamrTcPpgAEDiI6OBi70og4cOJDMzExGjBgBQEREBJmZmTgcDk6fPk1AQACenp4EBweTl5fnsqxvz5eRkQHAJ598wrBhw4ALd3X65JNPOHz4MIMHD3bWcPTo0VZZXxERERG5PLc453TVqlXs2bOH++67D5vNhp+fHwC+vr5YrVasVismk8nZ3tfXF5vN5rIMm82Gr69vnefPnTtXZ75L215UXV3dYH1r1qwhLCzM5UdEREREml+bnnN60RNPPEFOTg6zZ8/GbDZTWVkJgN1uJzAwkMDAQKqqqpzt7XY7ZrPZZRlmsxm73V7neR8fnzrzmUwmZ9uLunbt2mB98fHxxMfHu0xTQBURERFpfm3ac3rw4EGKiooACA4O5uTJk1gsFnJycgDIysrCYrHg7e1NUFAQZ86coba2ltLSUoYOHeqyrEvny87OZty4cQCMGTOG48ePA5Cbm8uYMWOIjIwkPz8fgMLCQkaOHNkq6ysiIiIil+dhGIbRVi9++PBhtmzZwuDBg/niiy/44Q9/SHR0NElJSfTo0YPy8nISExOdo/U3b96Mj48PFouFmJgYKioqmDZtGmlpaRiGwbJlyzCZTJw6dYr58+djMpk4efIkqamp9OjRg+DgYJfR+nl5edhsNuLi4q54tH5YWBi5ubktsFVEREQuz8PDgzb8+hZpUW0aTr+rjRs30rVrV+67775Wf22FUxERaSsKp9KRucU5p1dr0qRJ9O3bt63LEBEREZFm4haj9a+WgqmIiIhIx9Kuw6mIiIiIdCwKpyIiIiLiNhRORURERMRtKJyKiIiIiNtQOBURERERt6FwKiIiIiJuQ+FURERERNyGwqmIiIiIuA2FUxERERFxGwqnIiIiIuI2FE5FRERExG0onIqIiIiI21A4FRERERG3oXAqIiIiIm5D4VRERERE3EaXtnzx3NxcNmzYwJAhQ8jJyWHOnDlkZmayZcsWfHx8ANi4cSMA+/fv54MPPsDhcHD33XczevRol2U5HA6Sk5MJCgrCarWSkJCAl5cXOTk5bN26lW7dujF27FhiYmIAeOWVV3A4HJSWljJr1ix69uzZuisvIiIiInV4GIZhtNWLZ2Rk0L17d8LDw9m9ezcZGRmMGjWKqKgo+vXr52xXW1vL1KlT2bZtG+fOneOhhx5i27ZtLsvaunUr1dXVxMbGsnLlSsLDw7ntttt44IEHSE1Nxd/fn6lTp7JlyxaKiopYuXIlq1evJj09nffee48FCxZcUe1hYWHk5uY2y3YQERG5Eh4eHrTh17dIi2rTw/rjxo0jPDwcuBBATSYTAJs2beKll17i7bffBuDLL7/kmmuuwcvLC19fXxwOBzabzWVZmZmZjBgxAoCIiAgyMzNxOBycPn2agIAAPD09CQ4OJi8vj4MHDzrbhoeHk5mZedk616xZQ1hYmMuPiIiIiDS/Nj2sf5FhGOzevZvf//73OBwOJk6ciNls5oknnqBPnz54enri6+vrbO/r64vNZsNsNjun2Ww2/Pz8nM9brVasVqsz8F46n81mcy7v4rTLiY+PJz4+3mWaAqqIiIhI83OLAVFr165l2rRpXHvttfTv398ZOiMjIzl8+DBmsxm73e5sb7fbXYIpgNlsprKy0vl8YGAggYGBVFVV1Znv0uXVtywRkY7Kw8MDDw+Pti5DRKRBjYbTjz/+mK+++oqioiL+8Ic/kJWV1awFvP7660RERDB27Fjef/99VqxYQU1NDQDFxcWEhIQwcOBAysvLqampwW634+3tXSdQWiwWcnJyAMjKysJiseDt7U1QUBBnzpyhtraW0tJShg4dSlRUlLNtdnY2FoulWddJRERERK5OowOiFi5cyJw5c5g/fz7jx4/n888/57nnnmuWF09PT2f27NkMHToUgO7duxMVFUV+fj69e/fm66+/5qmnnsLT05P9+/ezb98+HA4HU6ZM4YYbbuDYsWOkpKSQkpKCw+EgKSmJHj16UF5eTmJionO0/ubNm/Hx8cFisbiM1q+qquLUqVPMmTPnikfra0CUiLRHF3tNNZimfdOAKOnIGg2nL730Eg888AAzZ87kT3/6E6+88gpxcXGtVd9lJSUlMW7cOMaPH9/qr61wKiLtkcJpx6BwKh1ZowOivvzySx555BGmTp1KXl6e83C4O3jwwQfp27dvW5chIiIiIs2k0Z7Tc+fOkZ+fz4gRIygqKqK6uppBgwa1Vn1uSz2nItIeqee0Y1DPqXRkDQ6IevPNNwHw8fFxXhO0X79+fPbZZ61TmYiIiIh0Og32nI4dOxZ/f3+XaYZhYLfb+eijj1qlOHemnlMRaY/Uc9oxqOdUOrIGzzmdN28emZmZTJ482dlzahgGf/7zn1utOBERERHpXC57zmlNTQ3vvvsux44d4/bbb2f06NGtWZtbU8+piLRH6jntGNRzKh1ZowOi4MIfseXLl5Odnc369etboy63p3AqIu2RwmnHoHAqHVmjl5I6fPgwa9eupaioiBkzZrRGTSIiIiLSSTUYTg8ePMi6desoLy9nxowZ3HHHHXh4eHD8+HGGDRvWmjWKiIiISCfR4GH94cOHExERwa233uoy/dChQ2zYsKEVSnNvOqwvUpcOGbs/7aOOQYf1pSNrsOf0kUceITY2ts70ysrKFi1IRERERDqvBntOHQ4H3t7eV7Sw2tpaPvzwQ7744gscDgfBwcGMHz+ea665plmKdSfqORWpS71y7k/7qGNQz6l0ZA3eIerIkSPs2bMHAKvVyvz585k7dy6FhYX1tj927BgPPvggb731Fv/973+prKzk008/5dFHH2Xz5s0tU72IiIiIdCgNHtZ/+eWXmTZtGgArVqzA09OTO+64g+XLl7Ny5UqXtjU1Nezbt4/169fX29v67rvv8sUXX3Ddddc1c/kiIm1DPVciIi2jwXA6atQooqOjqampYc+ePezevZuAgAA+//zzOm29vLyYOXNmgy/y4x//uHmqFREREZEOrcHD+hcdPHiQoUOHEhAQAECXLo1eGtXFN998c3WViYiIiEin02DSPHv2LElJSezdu5eEhAQA3nrrLfLy8q7oBTZs2NBgr2pubi4bNmxgyJAh5OTkMGfOHHr16kVycjJBQUFYrVYSEhLw8vIiJyeHrVu30q1bN8aOHUtMTIzLsi7excrf35+SkhISExMxmUwUFxfz4osv0rNnT3r16sX9998PwM6dOykoKKCsrIzp06cTGhp6ReslIiIiIs2vwdH6NTU1fPDBB/Ts2ZNRo0YBF3pRg4ODGTBgQJ32EydOpLa21mWaYRjY7XY++uijel88IyOD7t27Ex4ezu7du8nIyGDEiBFUV1cTGxvLypUrCQ8P57bbbuOBBx4gNTUVf39/pk6dypYtW1zOb92/fz8HDhwgMTGRbdu2UVlZycMPP8zcuXOZMWMGw4YNY9q0aSxfvhwfHx9mzJjBpk2bKCwsZPHixaSmpl7RhtNofZG6OtNI8PZ6zmln2kcdWXt9/4k0RYM9p15eXkyYMMH5+JtvvmHUqFF069at3vY/+clPmDRpEkFBQc5phmHw5z//ucEXHzdunPP32tpaTCYTmZmZPPDAAwBERESQkZHB+PHjOX36tPPUguDgYPLy8ggPD3fOn5mZyYgRI5zzrVq1iocffphPPvnEeUersLAwPvnkE0wmE4MHDwZgwIABHD16tMEaRURERKT1NHjO6YoVK1i0aBFwIfhFR0czceJEduzYUW/7Rx99FLvdTkhIiPOnX79+TJ06tdEiDMNg9+7dPPzww9hsNvz8/ADw9fXFarVitVoxmUzO9r6+vthsNpdl2Gw2fH196zx/7ty5OvNd2vai6urqButbs2YNYWFhLj8iIiIi0vwaDKenTp1i/vz5AKSkpLBkyRI++OADMjMz621vMpm46aab6kwfMmRIo0WsXbuWadOmce2112I2m513obLb7QQGBhIYGEhVVZWzvd1ux2w2uyzDbDZjt9vrPO/j41NnvkvbXtS1a9cG64uPjyc3N9flR0RERESaX4PhdODAgXh5eWGz2cjPz2fSpEl06dKF/v37N7iwTz/9lFWrVvG///u/rFq1ik8//bTRAl5//XUiIiIYO3Ys77//PhaLhZycHACysrKwWCx4e3sTFBTEmTNnqK2tpbS0lKFDh7os59L5srOznacMjBkzhuPHjwMXBmCNGTOGyMhI8vPzASgsLGTkyJGN1ikiIiIiLa/Bc04vHhbfuXMnEyZMcJ5Ef+bMmXrbr1u3juPHjzN69Gj69u1LZWUlGzZsID09nccee6zeedLT01m3bh1Dhw7l1VdfpXv37qxcuZKkpCTWrl1LVVUVEydOBGDBggW88MIL+Pj48Nhjj+Ht7U1FRQXTpk0jLS2Nm2++mYyMDNauXevS6ztnzhxSU1Pp0aMHP/zhD53nxN57770sXboUm81GYmLiVW4+EREREWlODY7W/9vf/sby5cvx9vbmpZdeom/fvixYsIBu3bqxYMGCOu3Xrl3Lr3/96zrT16xZQ3x8fPNXDmzcuJGuXbty3333tcjyL0ej9UXq6kwjwdvraOnOtI86svb6/hNpigZ7Tu+44w7uuOMOl2kXB0jV59SpU5w9exZ/f3/ntLNnz/Kf//ynGcqs36RJk+jbt2+LLV9EREREWteV3e7pMqZOncp9992Hh4cHfn5+nD17Fg8PD55//vnmeok6FExFREREOpYGD+tfrcLCQsrKyujRo0e9F+vvKHRYX6SuznTIuL0eVu1M+6gja6/vP5GmaHC0/vz586mpqWlwAFRDBgwYwOjRo53B9OKoeBERERGRxjQYTgcNGoSXlxevv/66y/TNmzc3eeElJSXMnTv36qsTERERkU6lwcP68+fPp6ioiOLiYkJCQoALh4FKSkrYs2dPqxbpjnRYX6SuznTIuL0eVu1M+6gja6/vP5GmaHBA1OLFizl16hSvv/46999/P3Dhj9m2bdsaXFhJSQkHDhxwnnN6880306dPn+avWkREREQ6pCYNiLJarRQWFjJw4MA6tw29aNu2bWzfvp2IiAj8/f05e/Ys2dnZ/M///A/33HNPsxfe1tRzKlJXZ+qVa689V51pH3Vk7fX9J9IUjV5K6i9/+QsrVqygZ8+elJWVMWvWLCZPnlyn3YkTJ9iyZUud6cnJyc1TqYiIiIh0eI2G088++4y///3veHh4UFtby/PPP19vOLXb7XWmGYbBN9980zyVioiIiEiH12g4DQgIcB4G8vT0JCAgoN520dHRTJ48mQEDBjgvwl9UVERCQkLzViwiIiIiHVaj4fTcuXPMnj2bkJAQioqKCA0NrbfdrbfeSlRUFJ999plzQFRkZCS+vr7NXbOIiIiIdFBNGhD1z3/+k7y8PMLCwvjBD37QGnW5PQ2Iks6oscE0nWmwTXsdkNKZ9lFH1l7ffyJN0eBF+C81fvx4fvnLXzYpmMbFxQHwi1/84rtVJiIiIiKdTpPC6dXQf3QiIiIicqVaLJyKiIiIiFypRgdExcTEkJKSwvDhw1ukgCNHjvDUU08xa9YsJkyYwI4dO9iyZQs+Pj4AbNy4EYD9+/fzwQcf4HA4uPvuuxk9erTLchwOB8nJyQQFBWG1WklISMDLy4ucnBy2bt1Kt27dGDt2LDExMQC88sorOBwOSktLmTVrFj179myR9RMRERGRpms0nE6cONElmB4/fpxhw4Y1uuCLJ9035uTJk4wYMcJl2vLly+nXr5/zcW1tLcuXL2fbtm2cO3eOhx56qM5tVNPS0ggNDSU2NpaVK1eyd+9ebrvtNhYtWkRqair+/v5MnTqVW265haKiIrKysli9ejXp6emkpqayYMGCJtUrIiIiIi2n0cP6FRUVLF++nLS0NHbu3ElKSsrlF+h5YZFNDac/+tGP6kzbtGkTL730Em+//TYAX375Jddccw1eXl74+vricDiw2Wwu82RmZjpDbkREBJmZmTgcDk6fPk1AQACenp4EBweTl5fHwYMHnW3Dw8PJzMxsUq0iIiIi0rIa7TnNyckhJCSEkydPAlBeXn7Z9i+99BIAr776KgAlJSX06dOnyQXdeOONTJw4EbPZzBNPPEGfPn3w9PR0uV6qr68vNpsNs9nsnGaz2fDz83M+b7VasVqtmEymOvPZbDbn8i5Ou5w1a9Y0GspFRERE5LtrNJwuWrSIiIgI5+MTJ040eeG1tbWkpKTwhz/8ocnz9O/f3/l7ZGQkhw8fZuLEiS63R7Xb7S7BFMBsNlNZWel8PjAwkMDAQKqqqurMZzabsVqtDS7r2+Lj44mPj3eZFhYW1uR1EhEREZGmafSwfmBgIHPnzmXhwoX84x//aLCXcceOHVgsFiZPnkxxcTHbt2+v95B9Y1asWEFNTQ0AxcXFhISEMHDgQMrLy6mpqcFut+Pt7V0nUFosFnJycgDIysrCYrHg7e1NUFAQZ86coba2ltLSUoYOHUpUVJSzbXZ2NhaL5YrrFBEREZHm12jP6YsvvsjPf/5zDh48yPjx41m8eDHXX399nXb//Oc/+cc//kFZWRlz585lyJAhrF+/nr59+152+W+99Ra5ubnU1NQQFBREYGAgzzzzDL1796ampoY77rgDT09P5syZw+LFi3E4HDz55JMAHDt2jJSUFFJSUpgyZQpJSUmsXbuWqqoqJk6cCMCCBQt44YUX8PHx4bHHHsPb25vBgwdz/fXXs3r1ak6dOsWcOXOuZtuJiIiISDNrNJwOGjSIqKgoPvvsM7y9vbn22mvrbTds2DACAgIICAhg9OjRJCYmAnDq1KkG5wG48847ufPOO52P6wu+ANHR0URHR7tM27VrF/fccw8A3t7ePPXUU3XmGzFiBM8++2yd6RfvZCUiIiIi7qPRcJqXl8d7771HWVkZ6enpzoFR33by5EnS09MBOHv2LBkZGRiGwdtvv83zzz/fvFX//x588MFGe2ZFREREpP3wMBq5z2hpaSlJSUnk5eURFhZGQkJCvT2hkyZNqjcolpSUsGfPnuar2E2EhYWRm5vb1mWItKqLl4hr6M9GY893JB4eHu1yPTvTPurI2uv7T6QpGg2ncGHU/ddff03Pnj2d1zH9tvfee49JkyY1eXp7p3AqnZHC6f9pr+GgM+2jjqy9vv9EmqLRcLpnzx6efvppPD09qa2t5emnn3beArQ+R44coaioiH79+jFq1KhmL9hdKJxKZ9SZwmljX/7tNRx0pH3UmbXX959IUzR6zumuXbvYvXs3/v7+VFRUkJiYWG84PXPmDL/+9a+prq6mT58+lJSU4OXlxbp167jmmmtapHgRERER6VgaDaeRkZH4+/sD0L17d0aOHFlvu9TUVH7zm99w4403OqcdPHiQ1NRUfve73zVTuSLuQb0WIiIiLaPBcHrxdp3Z2dn8+9//dt7C9Ouvv663vclkcgmmAFFRUc4R/CIiIiIijWkwnObk5BATE0NISIhzWkhICO+//3697b29veud7uPj8x1LFBEREZHOosEBUf/5z3/o3bt3nemffvopN9xwQ53pEydOpE+fPi7TDMOgtLSUvXv3NlO57kMDojq3znpYXwOimv68u+pI+6gza6/vP5GmaLDn9GIwPXnyJP/4xz84e/YsAIcOHWLDhg112t90003cdddddabv2rWrmUoVERERkY6u0QFRc+fOZfLkyc7D+8ePH6+33e9+9zu6d+9eZ/rw4cO/Y4kiIiIi0lk0Gk5HjhzJtGnTnI8bunbpE088wS9+8QtGjx4NgK+vL56engQEBDRTqSIiIiLS0TUaTm+88Ubmz5/v7Dlt6LD+jTfeSHR0NCkpKezfv58XXniBfv36NXvBIiIiItJx1X8v0kts2LCBESNGEBISQkhISIM9oRdPsn/88ccZNWqUM5ieOnWqGcsVERERkY6s0Z7TG264oUmH9cvLyzlx4gQAZ8+edf6+adMmEhMTm6NWkTajkbEiIiKto8FLSV3029/+li5dujh7Qhs6rB8VFUX37t3rfIFXVlaSmZnZfBW7CV1KqnP5djjtrGFVl5Jq+vPuqiPto86svb7/RJqi0Z7T4uJi7rnnHufj+kbrG4bBk08+yd13313nuZ07dwJQVVWFyWSq8/yRI0d46qmnmDVrFhMmTMDhcJCcnExQUBBWq5WEhAS8vLzIyclh69atdOvWjbFjxxITE1OnhuXLl+Pv709JSQmJiYmYTCaKi4t58cUX6dmzJ7169eL+++931lVQUEBZWRnTp08nNDS0sU0hIiIiIi2s0XNOly5dypQpU5w/Tz/9dJ02Hh4eZGdn19uTOHHiRBYvXozVaq13+SdPnmTEiBHOx2lpaYSGhjJz5kxMJpPzAv6LFi1izpw5JCQksG7dOhwOh8tyDhw4wPnz55kxYwbh4eFs3boVgGXLlhEbG8sTTzzB7t27OX36NBUVFbzxxhvMnj2buLg4kpKSGtsMIiIiItIKGu05/eijj/joo4+cj/fu3cvq1avrtJs1axaLFi3io48+wmw207VrV2w2G9dccw0LFy6kb9++9S7/Rz/6Ef/85z+djzMzM3nggQcAiIiIICMjg/Hjx3P69GnnYKzg4GDy8vIIDw93me9iyI2IiGDVqlU8/PDDfPLJJwwbNgy4cCj+k08+wWQyMXjwYAAGDBjA0aNHG9sMIiIiItIKGg2n27dvx2KxAFBaWkqXLvXP4u/vz5IlS6isrOTEiRM4HA6Cg4PrvQXq5dhsNvz8/IAL10q1Wq1YrVaXUwJ8fX2x2Wx15vP19a3z/Llz5+rMd+7cOWfbi6qrq+natWu9Na1Zs4aUlJQrWg8R0XlxIiJy5RoNp0uWLHFe4xTgz3/+82Xb+/n5fae7QpnNZiorKwGw2+0EBgYSGBhIVVWVs43dbsdsNteZz26313nex8enznwmk8nZ9qKGgilAfHw88fHxLtPCwsKuYu1ERERE5HIaPee0uLjYeWh/37597N+/v0ULslgs5OTkAJCVlYXFYsHb25ugoCDOnDlDbW0tpaWlDB06tMH5srOzGTduHABjxoxxDuLKzc1lzJgxREZGkp+fD0BhYSEjR45s0XUSERERkaZp9FJSd999t/NcTj8/P376058SGRnZbAW89dZbrF+/nuuuu45p06YRFhZGUlISPXr0oLy8nMTEROdo/c2bN+Pj44PFYiEmJoaKigqmTZtGWloahmGwbNkyTCYTp06dYv78+ZhMJk6ePElqaio9evQgODjYZbR+Xl4eNpuNuLi4Kx6tr0tJdS66lNQFV3opqfa8nXQpKXFn7fX9J9IUjYbT7OxsIiIiWqueK7Jx40a6du3Kfffd1+qvrXDauSicXqBw2vTn3ZXCacfQXt9/Ik3R4GH99evXA7htMAWYNGlSmwRTEREREWkZDQ6Ievvtt+tccP/DDz/EMAw++OCDFi+sKRq6PJVIW1OvhoiIyNVpMJzOmzePm266CbhwOaZFixYxaNAgli1b1mrFiYi4C/3D0f7oFAaR9qnBw/oXg2l+fj4///nP6d27N+vXr6dnz56tVpyIiIiIdC6Xvc7pjh07SElJ4dlnnyU6Orq1ahIRN6GeJxERaW0N9pwmJiby5ptvsnnzZpdgumHDhtaoS6RduRjiRNyB3o8i0p412HOal5dHWFgYK1eudJmem5vLww8/3NJ1iYiIiEgn1GA4/e1vf8v3v//9OtM//PDDFi1IRETcl071EJGW1uhF+KV+ugh/59LYRfg76kX667uo/qWPm9K+vW6H9ryPL1fbdw2X7Smctqdar5Q7v/9EvqsGzzkVEREREWltCqciIiIi4jYUTkVERETEbSicirQwDw8PXdpHRESkiRRORURERMRtKJyKiIiIiNtQOBURERERt9HgRfjb2sSJEwkJCQEgLi6OyMhIli9fTnBwMOfPnyc+Pr7OPDt37qSgoICysjKmT59OaGgoDoeD5ORkgoKCsFqtJCQk4OXlRU5ODlu3bqVbt26MHTuWmJiY1l5FEREREfkWt+05nTJlChs3bmTjxo3ccsstvPrqq9x66608+uijFBYWcvToUZf2FRUVvPHGG8yePZu4uDiSkpIASEtLIzQ0lJkzZ2Iymdi7dy8AixYtYs6cOSQkJLBu3TocDkerr2NTaDCNiIiIdCZuG04PHTrEH//4R9atW8fZs2fJzMxk5MiRAISHh5OZmenS/vDhwwwePBiAAQMGOMNrZmYmI0aMACAiIoIzNYm+AAAR4ElEQVTMzEwcDgenT58mICAAT09PgoODycvLa8W1ExEREZH6uO1h/Xnz5hEREcG+fftYsmQJNpsNX19fAPz8/CgqKnJpf+nzF1VXV2Oz2fDz8wPA19cXq9WK1WrFZDI52/n6+mKz2RqsZc2aNaSkpDTXqom4rY58u0cREWkf3LbnNCIiAoDIyEg+/fRTzGYzdrsdgMrKSsxms0v7S5+/qGvXrpjNZiorKwGw2+0EBgYSGBhIVVWVs53dbq+zvEvFx8eTm5vr8iPSUei0ERERcSduGU7T09M5cOAAAMXFxYSEhGCxWJyH6rOzs4mKinKZJzIykvz8fAAKCwudpwBYLBZycnIAyMrKwmKx4O3tTVBQEGfOnKG2tpbS0lKGDh3aWqsnIiIiIg3wMNzw+N2xY8dISUkhIiKCgoIC4uLi6NWrF8uWLaN3794YhuEcrT9lyhRee+01AgIC2LlzJ3l5edhsNuLi4pyj9ZOSkujRowfl5eUkJiY6R+tv3rwZHx8fLBbLFY/WDwsLa5UeVB1mdQ8eHh4u++BKHrv7PrxcrY09rm9Z327vruvdmCvd5+7kcrV91/eju7+fL9Weam1MR/psiTTGLcNpU33++ee89tprJCcnt/prK5x2LgqnCqcKp80zf2tqT7U2piN9tkQa45aH9ZuqV69e/OEPf2jrMlqVLi0lIlI//X0U6RjcdrR+U/Tu3butSxARERGRZtSue05FRKTzUk+pSMekcCoiIiIibkPhVERERETchsKpiIiIiLgNhVMRERERcRsKpyLSKXWmwTSdZT1FmlNn+hvhbhRORURERMRtKJyKiIiIiNtQOJVWo8MjoveASMvQZ0s6EoXTFqA/EiIiIiJXR+FURKQZ6J9S96cBLiLtg8KpiLQZBYW2p8AmIu5G4dTNdOYvis663iIi0vF15u/3K6VwKiJXrSP9oW1PXxwtWac7b4e2rq05X9tdt7E7uHTbtPU+l6vzXfdZl2aqo13auXMnBQUFlJWVMX36dEJDQ9u6pE7Fw8MDwzBarP13cfGD1Vqv1xF9+49Ta3/BaB82D21HaWlX8re9Pb0f21Ot31Vzr2un7TmtqKjgjTfeYPbs2cTFxZGUlNQmdXzX/wpb8wu/sVq//Xx7WrfOpLm3a2PL64z7sbV7ezryNu7I6yYdy5V+7ttzr3BL195pw+nhw4cZPHgwAAMGDODo0aMt9lo6ROH+WnO/NBbim/K4seVf7nFT5mkPOvIXwZXuc3fyXf8hbctQ39yH7a/kb39L/3P/Xda1pWtpT9pz7d9Fq382jc7Q31yPt99+m6ysLJ588kkAxo8fz549e+jatWudtmvWrCElJaW1SxQRERFpt3Jzc69qvk57zqnZbMZut7tMqy+YAsTHxxMfH+8yLSws7Ko3ukhHos+CyAX6LIj8n7CwsKuet9Me1o+MjCQ/Px+AwsJCRo4c2cYViYiIiEin7TkNCAjg3nvvZenSpdhsNhITE9u6JBEREZFOr9OGU4C77767rUsQERERkUt4Pf3000+3dRHtlcViaesSRNyCPgsiF+izIPJ/rvbz0GlH64uIiIiI++m0A6JERERExP0onIqIiIiI21A4FRERERG3oXAqIiIiIm6jU19K6mosXLiQgoICAIYPH87vf/97AF555RUcDgelpaXMmjWLnj17tmWZIq1i586dFBQUUFZWxvTp0wkNDW3rkkRazcSJEwkJCQEgLi6OyMhIli9fTnBwMOfPn69zZ0GRjuTIkSM89dRTzJo1iwkTJuBwOEhOTiYoKAir1UpCQgJeXl7k5OSwdetWunXrxtixY4mJiWl02QqnV6hXr148++yzLtPy8/PJyspi9erVpKenk5qayoIFC9qoQpHWUVFRwRtvvMGmTZsoLCxk8eLFpKamtnVZIq1mypQpLgF0+fLl3HrrrUyYMIF58+Zx9OhR3X1QOqyTJ08yYsQI5+O0tDRCQ0OJjY1l5cqV7N27l9tuu41FixaRmpqKv78/U6dO5ZZbbsHb2/uyy9Zh/StUWVnJyy+/zKpVqzhy5AgABw8edO6g8PBwMjMz27JEkVZx+PBhBg8eDMCAAQM4evRoG1ck0roOHTrEH//4R9atW8fZs2fJzMx0hlF9F0hH96Mf/cjlcWZmpjMLRUREkJmZicPh4PTp0wQEBODp6UlwcDB5eXmNLls9p1fozjvvZPjw4VRXV3PPPfewadMmbDYbvr6+APj6+mKz2dq4SpGWd+n7/qLq6mq6du3aRhWJtK558+YRERHBvn37WLJkictnws/Pj6KiojauUKT12Gw2/Pz8gAtZyGq1YrVaMZlMzjZNzUgKp/X473//y29+85s603/1q18xYcIEAHx8fOjbty/5+fmYzWasVisAdrsds9ncqvWKtAWz2YzdbneZpmAqnUlERAQAkZGRLF261PmZ6N69O5WVlfoukE7FbDZTWVkJXMhCgYGBBAYGUlVV5WzT1IykcFqPXr16sXnz5nqfW7p0KfPmzQPg1KlTBAcH4+/vz8qVKwHIzs7W7eukU4iMjHSeY1pYWKhz66RTSU9Pp7a2lptvvpni4mJCQkIYPnw4R48e5dprryU7O5tHHnmkrcsUaTUWi4WcnBy+973vkZWVhcViwdvbm6CgIM6cOYO/vz+lpaUMHTq00WXp9qVXaN68efTt25eqqioGDhzIAw88AFwYrV9VVcWpU6eYM2eORutLp7Bz507y8vKw2WzExcVptL50GseOHSMlJYWIiAgKCgqIi4ujV69eLFu2jN69e2MYhkbrS4f21ltvsX79eq677jqmTZtGWFgYSUlJ9OjRg/LychITE52j9Tdv3oyPjw8Wi6VJo/UVTkVERETEbWi0voiIiIi4DYVTEREREXEbCqciIiIi4jYUTkVERETEbSicioiIiIjbUDgVEZF24/z58606n4i0Pq+nn3766bYuQkTk2w4ePMjjjz/Ov/71LwoKCtiyZQuGYXDddde1dWnthsPh4Pe//z179uwhJiaGY8eOsWbNGiZMmIBhGMyfP58vvviCN998Ez8/P1JTU9m/fz9eXl4MHDiwrct3UV1dzdKlSxk+fDjp6en88pe/5NixY+Tm5vLOO+/g6el52evsnj59mpSUFCwWC15eXq1XuIhcMd0hSkTcUlRUFGFhYXzve9/jnnvuwWq1cvvtt/PDH/6wrUtrN7y9vZkyZQppaWkADB8+nGeeeQaA//znP5w4cYIlS5ZQXV3NM888w09/+lPGjRtHTU1NW5Zdr9TUVMaOHUvv3r2JiYnhtdde46c//Snf//73OX/+PDNmzMDLy4vo6Oh65+/duzc33ngjqamp9d6eWkTch8KpiLQLZWVlBAYGkpWVxYIFCxg5ciT+/v6899577N27lxUrVtClSxcqKysZOXIkd955J2fOnGHp0qUMGDCAwsJCbrzxRry9vVmwYAGHDh0iMzOT+fPn8/rrr5Ofn89zzz3H+PHjOX/+PDk5OcybN48dO3YwaNAgsrOzWbVqFeXl5SQnJxMaGspXX31FbGws/v7+zJs3j6CgIPr168c777zDu+++S0BAAAC7du0iOTmZhx56iLy8PGpqarjpppv48MMP6dq1K8nJyQA8+uijjB49mi+//JIf//jHREdHc+jQoTo17N69m4yMDPr06cPJkyd59tlnXbZVfn4+ycnJjBo1irKyMuf05ORksrKy2LhxI6+99holJSWsWbOGm2++maysLM6fP09FRQXjxo1jyZIlDBo0iMLCQv7f//t/BAQE1LuOu3btoqSkBG9vb7p378706dNZsWIFu3btIjY2lkOHDhEWFsbs2bMBWLlyJV5eXpw/f55vvvmG+fPnk5mZyTvvvEP//v0pLi7mySefxNvb22Wdtm/fTlxcXL3vjS5duvDII4+wfv16oqOjWbZsGZ6enjgcDgIDA/nVr34FwM0338wzzzyjcCri7gwRETeVmJhoPP7440Zqaqrx/PPPG3l5eYZhGMbq1auNpUuXGoZhGJ9//rmxb98+49FHHzUMwzBqa2uNyZMnG+Xl5UZycrLxpz/9yTAMw6iqqjLeffddwzAMY8KECc7XiI2NNU6cOOF8vU2bNjmX+9xzzxmpqalGbW2t8fHHHxuGYRhJSUnOZf773/82pkyZYhiGYWzfvt2YNWuWYRiGkZubazgcDpd1iY2NNT788EPDMAzjrrvuMtLT052/l5WVGYZhGHv37jUMwzBsNptzufXVMHPmTGPXrl2GYRjOaZd69NFHjffff98wDMN44403jMTERMMwDOPEiRNGbGxsnd8vrntGRoZhGIaRnJxsrF+/3tnu5z//eb3rmJOTY9x9990u6/jvf//bMAzDuP76642Kigrj/PnzRnR0tGEYhrFv3z7jN7/5jbP99u3bjdraWuMHP/iBcxusXr3a2LJli8v6nDlzxoiKiqqzPQ8cOOB8/MUXXxi33367y3a8uH0rKiqcj8eNG2dYrdY620xE3Id6TkXErd1yyy3cc889daYPGjQIgJEjR/LHP/6RiooKXn75ZQCGDBlCWVkZx48fJyoqCoBu3boxefLkRl/v0uX27t2bdevW8bOf/YxbbrmFG264gby8PCwWCwD9+vXj+PHjdeYdNmxYvcvu378/AN27d6dfv37O3ysrK+nevTvHjh3j888/x9vbG5vNBsDMmTPr1PC73/2OF198kQ0bNnDXXXcxZswYl9cpKChwvla/fv34+OOPG13vS+Xl5VFWVsbLL7+MYRiYzeZ61/Gvf/0rNTU1zu0eHBxMWVkZgwcPpmfPnvj7+wPQtWtXAI4fP05ISIhzWT/72c8oKyujoqKCbdu2AVBeXk5gYKBLPQ6Hgy5dLv91VVxc7NymdrudF154gYCAACoqKrDZbM5aunTpwrlz565oe4hI61I4FZF2ycPDw/l7WFgYOTk5zsO3f//73wkODiYsLIyioiIAzp49ywcffMDkyZPx8/Pj7Nmz+Pv7U1pa2uBys7KyWLhwIefPn+fBBx/k6NGjLss8ceKESxC9dN4rtW/fPj799FNefvllqqur2bJlS4M1nDp1isWLF1NVVcVPfvIT7rrrrjoB8quvvmLIkCHOWq9EWFgY/fv3595778UwDN55551613HYsGH4+/s7t3t6erpzIFV92yIsLIw333zT+Xjbtm1MnTqVwMBAYmNj8fX15auvvnIG84sCAwM5d+4chmHUu9zz58+zYcMGHnroIc6cOcPChQs5ePAgXl5e7N2719nOMAzOnTtHUFDQFW8TEWk9Gq0vIm7p4MGD7Ny5k6+//pqQkBCCg4OBC4Fw48aNlJSUMHjwYIKCghg4cCB5eXkcOHCAjIwMqqqqsFgsXH/99aSlpZGXl8fevXuZMGECgYGBGIbB9u3b+e9//0tubi52u52AgAC2bdvG119/zfDhwwkICOBf//oXf/nLXzh69CheXl787Gc/Y/To0ezYsYPjx4/z/vvvM3fuXLy8vFi/fj35+fn069fPpXcQICMjg507d9KlSxc8PDzYsWOHy+8At912G++88w65ubnk5uaSnp7O4MGDOXnyZJ0a3nrrLTIyMvjss8/o06cPt99+u8vrDR8+nNTUVAoKCvjyyy85evQow4cP55133uHjjz/muuuu4+9//zsff/wxXbp0oba2lu3bt/P111/Tv39/JkyYQFpaGseOHeNvf/sbQ4cO5Zprrqmzjj169OCbb77hr3/9K4cOHeLEiRNMmjSJtLQ03nvvPa677joKCgp499136d27N3fccQfHjh0jIyODjIwMBgwYwJAhQwgPD+fVV1/l888/58CBA9x+++2YTCbn+nh4eJCdne3c33v27OEvf/kLlZWV5Obm8u6773LnnXdy66234u3tzfHjx9m7dy8nTpzgww8/xNPTk6ioKPLy8iguLubHP/5xC797ReS78DAMw2jrIkRERC6nqKiIlJQUnnvuOedpAleiurqahQsX8vjjj9f550FE3IvCqYiItAtlZWXU1tZe1WH506dP4+npSY8ePVqgMhFpTgqnIiIiIuI2dPtSEREREXEbCqciIiIi4jYUTkVERETEbSicioiIiIjbUDgVEREREbfx/wHKv0i7RGChrgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 756x207.664 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "width = 7\n",
    "height = width / 1.618\n",
    "fig, ax = plt.subplots(figsize=(width * 1.5, height / 1.5))\n",
    "\n",
    "# Exclude unmodified SSMs.\n",
    "mask = mass_groups['mass_diff_median'].abs() > tol_mass\n",
    "ax.bar(mass_groups[mask]['mass_diff_median'], mass_groups[mask]['num_ssms'],\n",
    "       width=0.4, log=False, color='black')\n",
    "\n",
    "ax.set_xlim((-50, 100))\n",
    "\n",
    "# Format y-axis numbers.\n",
    "ax.yaxis.set_major_formatter(mticker.StrMethodFormatter('{x:,g}'))\n",
    "\n",
    "sns.despine(ax=ax)\n",
    "    \n",
    "# Set tick labels at nice positions.\n",
    "ax.xaxis.set_major_locator(mticker.MultipleLocator(50))\n",
    "ax.set_xlabel('Precursor mass difference (Da)')\n",
    "ax.set_ylabel(f'Number of SSMs\\n(FDR={max_fdr})')\n",
    "\n",
    "plt.savefig('mass_diff.pdf', dpi=300, bbox_inches='tight')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
